/*============================================================================
  MakeLCP.cpp  -- Jun10, 2013 -- J.Benoit
  This program was derived from lcpcreate, which was originally bundled with
  the FARSITE 4 source code. It's purpose is to create Landscape (.LCP) files,
  the default format for FARSITE. This rewrite was intended to be more robust
  than the original. More input validation and error-checking is performed.
*/
#include<fstream>
#include<iostream>
#include<cmath>
#include<string>

using namespace std;

#define NODATA -9999

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ASCIIGridHeader struct
  Holds data from an ESRI Arc ASCII Grid layer file.
*/
typedef struct
{
  long   NRows, NCols;
  double XLLCorner, YLLCorner, CellSize;
} ASCIIGridHeader;

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  HeadData
  Header for landscape file.
*/
typedef struct
{
  long   CrownFuels;    // 20 if no crown fuels, 21 if crown fuels exist
  long   GroundFuels;   // 20 if no ground fuels, 21 if ground fuels exist
  long   Latitude;
  double loeast, hieast, lonorth, hinorth;
  long   loelev, hielev;
  long   numelev;       // -1 if more than 100 categories
  long   elevs[100];
  long   loslope, hislope;
  long   numslope;      // -1 if more than 100 categories
  long   slopes[100];
  long   loaspect, hiaspect;
  long   numaspect;     // -1 if more than 100 categories
  long   aspects[100];
  long   lofuel, hifuel;
  long   numfuel;       // -1 if more than 100 categories
  long   fuels[100];
  long   locover, hicover;
  long   numcover;      // -1 if more than 100 categories
  long   covers[100];
  long   loheight, hiheight;
  long   numheight;     // -1 if more than 100 categories
  long   heights[100];
  long   lobase, hibase;
  long   numbase;       // -1 if more than 100 categories
  long   bases[100];
  long   lodensity, hidensity;
  long   numdensity;    // -1 if more than 100 categories
  long   densities[100];
  long   loduff, hiduff;
  long   numduff;       // -1 if more than 100 categories
  long   duffs[100];
  long   lowoody, hiwoody;
  long   numwoody;      // -1 if more than 100 categories
  long   woodies[100];
  long   numeast, numnorth;
  double EastUtm, WestUtm, NorthUtm, SouthUtm;
  long   GridUnits;     // 0 for metric, 1 for English
  double XResol, YResol;
  short  EUnits, SUnits, AUnits, FOptions, CUnits, HUnits, BUnits, PUnits,
         DUnits, WOptions;
  char   ElevFile[256], SlopeFile[256], AspectFile[256], FuelFile[256],
         CoverFile[256], HeightFile[256], BaseFile[256], DensityFile[256],
         DuffFile[256], WoodyFile[256], Description[512];
} HeadData;

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  CellData
  Structure for holding basic cell information.
*/
typedef struct
{
  short e;                 // elevation
  short s;                 // slope
  short a;                 // aspect
  short f;                 // fuel models
  short c;                 // canopy cover
} CellData;

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  CrownData
  Structure for holding optional crown fuel information.
*/
typedef struct
{
  short h;  // canopy height
  short b;  // crown base
  short p;  // bulk density
} CrownData;

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  GroundData
*/
typedef struct
{
  short d;  // duff model
  short w;  // coarse woody model
} GroundData;

bool   GenerateLCP( string& ElevFileName, string& SlopeFileName,
                    string& AspectFileName, string& FuelFileName,
                    string& CoverFileName, string& TreeHeightFileName,
                    string& CrownBaseHeightFileName,
                    string& CrownBulkDensityFileName,
                    string& DuffLoadingFileName, string& WoodyFileName,
                    long Latitude, string& LandFileName );
void   SkipASCIIHeader( fstream& ASCIIFile, ASCIIGridHeader& Header,
                        const string& FileName1, const string& FileName2 );
void   FillCats( CellData *cell, CrownData *cfuel, GroundData *gfuel,
                 long AllCats[][100], long *NumAllCats,
                 double *maxval, double *minval );
void   SortCats( long AllCats[][100],long *NumAllCats, double *minval );
void   FinalizeHeader( HeadData *NewHeader, long AllCats[][100],
                       long *NumAllCats, double *maxval, double *minval,
                       string& ElevFileName, string& SlopeFileName,
                       string& AspectFileName, string& FuelFileName,
                       string& CoverFileName, string& TreeHeightFileName,
                       string& CrownBaseHeightFileName,
                       string& CrownBulkDensityFileName,
                       string& DuffLoadingFileName, string& WoodyFileName );
void   WriteHeader( fstream& OutFile, HeadData& Header );
void   WriteCell( fstream& OutFile, CellData& Cell );
void   WriteCrown( fstream& OutFile, CrownData& Crown );
void   WriteGround( fstream& OutFile, GroundData& Ground );

/*=============================================================================
  Main
*/
int main( int argc, char *argv[] )
{ //main
  long   Latitude = 35;     //Default latitude
  string LandFileName = "", ElevFileName = "", SlopeFileName = "",
         AspectFileName = "", FuelFileName = "", CoverFileName = "",
         OutputFileName = "";
  string TreeHeightFileName, CrownBaseHeightFileName, CrownBulkDensityFileName,
         DuffLoadingFileName, WoodyFileName;

  if( argc <= 1 ) {
    cout << "MakeLCP makes a binary FARSITE Landscape (.LCP) file.\n\n"
         << "All GIS grid files must be coregistered and with identical \n"
         << "resolution with the correct units. Preceed each file\n"
         << "with the identifying switch indicated below:\n\n"
         << "  Required inputs:\n"
         << "    SWITCH         NAME               (UNITS)\n"
         << "    -latitude      latitude           (-90 to 90)\n"
         << "    -landscape     Landscape                     (output file)\n"
         << "    -elevation     Elevation          (m)        (input file)\n"
         << "    -slope         Slope              (deg)      (input file)\n"
         << "    -aspect        Aspect             (deg)      (input file)\n"
         << "    -fuel          Fuel                          (input file)\n"
         << "    -cover         Cover              (%)        (input file)\n\n"
         << "  Optional GIS Themes:\n"
         << "    SWITCH         NAME               (UNITS)\n"
         << "    -height        Tree Height        (m)        (input File)\n"
         << "    -base          Crown Base Height  (m)        (input File)\n"
         << "    -density       Crown Bulk Density (kg/m3)    (input File)\n"
         << "    -duff          Duff Loading       (Mg/Ha)    (input File)\n"
         << "    -woody         Woody Profile                 (input File)\n"
         << endl;
    exit(0);
  }

  //Parse arguments.
  int i = 0;
  while( ++i < argc ) {
    if( argv[i][0] == '-' ) {  //Ensure keywords are lowercase
      for( int j = 1; j < strlen( argv[i] ); ++j )
        argv[i][j] = tolower( argv[i][j] );
    }

    if( ! strncmp("-lan", argv[i], 4) ) {
      if( argv[i+1][0] != '-' ) {
        LandFileName = argv[++i];

        //Make sure there isn't already a '.lcp' on file name.
        if( LandFileName.length() > 4 ) {
          string tail = LandFileName.substr( LandFileName.length()-4 );
          if( ! (tail[0] == '.' && (tail[1] == 'l' || tail[1] == 'L') &&
                 (tail[2] == 'c' || tail[2] == 'C') &&
                 (tail[3] == 'p' || tail[3] == 'P')) )
            LandFileName += ".lcp";
        }
        else LandFileName += ".lcp";
      }
      else {
        //User may have forgotten to follow switch with a file name.
        //This parsing loop will interpret the next arg as a switch.
        cerr << "## Odd file name following " << argv[i] << " switch: "
             << argv[i+1] << " ##" << endl;
      }
    }
    else if( ! strncmp("-ele", argv[i], 4) )
      ElevFileName = argv[++i];
    else if( ! strncmp("-slo", argv[i], 4) )
      SlopeFileName = argv[++i];
    else if( ! strncmp("-asp", argv[i], 4) )
      AspectFileName = argv[++i];
    else if( ! strncmp("-fue", argv[i], 4) )
      FuelFileName = argv[++i];
    else if( ! strncmp("-cov", argv[i], 4) )
      CoverFileName = argv[++i];
    else if( ! strncmp("-lat", argv[i], 4) )
      Latitude = atol( argv[++i] );
    else if( ! strncmp("-hei", argv[i], 4) )
      TreeHeightFileName = argv[++i];
    else if( ! strncmp("-bas", argv[i], 4) )
      CrownBaseHeightFileName = argv[++i];
    else if( ! strncmp("-den", argv[i], 4) )
      CrownBulkDensityFileName = argv[++i];
    else if( ! strncmp("-duf", argv[i], 4) )
      DuffLoadingFileName = argv[++i];
    else if( ! strncmp("-woo", argv[i], 4) )
      WoodyFileName = argv[++i];
  }

  while( LandFileName.length() == 0 ) {
    cout << "Output Landscape file: ";
    cin >> LandFileName;
  }
  while( ElevFileName.length() == 0 ) {
    cout << "Elevation file: ";
    cin >> ElevFileName;
  }
  while( SlopeFileName.length() == 0 ) {
    cout << "Slope file: ";
    cin >> SlopeFileName;
  }
  while( AspectFileName.length() == 0 ) {
    cout << "Aspect file: ";
    cin >> AspectFileName;
  }
  while( FuelFileName.length() == 0 ) {
    cout << "Fuel file: ";
    cin >> FuelFileName;
  }
  while( CoverFileName.length() == 0 ) {
    cout << "Cover file: ";
    cin >> CoverFileName;
  }

  cout << "Combining grid layers..." << endl;

  if( ! GenerateLCP(ElevFileName, SlopeFileName, AspectFileName,
                    FuelFileName, CoverFileName, TreeHeightFileName,
                    CrownBaseHeightFileName, CrownBulkDensityFileName,
                    DuffLoadingFileName, WoodyFileName,
                    Latitude, LandFileName) )
    cerr << "## MakeLCP Failed: Probable file read error - check paths ##\n";
  else
    cout << "Completed.\n"
         << "Landscape File Generated: " << LandFileName << endl;

  return 0;
} //main

/*=============================================================================
  GenerateLCP
  Generate the LCP file.
  Note: Units in this file are assumed, these must be determined elsewhere
  from user input.
*/
bool GenerateLCP( string& ElevFileName, string& SlopeFileName,
                  string& AspectFileName, string& FuelFileName,
                  string& CoverFileName, string& TreeHeightFileName,
                  string& CrownBaseHeightFileName,
                  string& CrownBulkDensityFileName,
                  string& DuffLoadingFileName, string& WoodyFileName,
                  long Latitude, string& LandFileName )
{ //GenerateLCP
  /*Open input files and read in & verify grid headers.
    Resolution and coordinate info is taken from Elevation file only.
  */
  ASCIIGridHeader ElevHeader;
  string keyword, value;
  fstream ElevFile( ElevFileName.c_str() );
  if( ! ElevFile.good() ) {
    cerr << "## GenerateLCP: Cannot open Elevation file " << ElevFileName
         << " ##" << endl;
    return false;
  }
  else {
    for( int i = 0; i < 6; ++i ) {
      ElevFile >> keyword >> value;
      if( keyword[0] == 'n' || keyword[0] == 'N' ) {
        if( keyword[1] == 'c' || keyword[1] == 'C' )
          ElevHeader.NCols = atoi( value.c_str() );
        else if( keyword[1] == 'r' || keyword[1] == 'R' )
          ElevHeader.NRows = atoi( value.c_str() );
      }
      else if( keyword[0] == 'x' || keyword[0] == 'X' )
        ElevHeader.XLLCorner = atof( value.c_str() );
      else if( keyword[0] == 'y' || keyword[0] == 'Y' )
        ElevHeader.YLLCorner = atof( value.c_str() );
      else if( keyword[0] == 'c' || keyword[0] == 'C' )
        ElevHeader.CellSize = atof( value.c_str() );
      else if( keyword[0] == 'n' || keyword[0] == 'N' ) {}
        //Skip if 'NO_DATA_VALUE' keyword.
    }
  }

  fstream SlopeFile( SlopeFileName.c_str() );
  if( ! SlopeFile.good() ) {
    cerr << "## GenerateLCP: Cannot open Slope file " << SlopeFileName
         << " ##" << endl;
    return false;
  }
  else
    SkipASCIIHeader( SlopeFile, ElevHeader, SlopeFileName, ElevFileName );

  fstream AspectFile( AspectFileName.c_str() );
  if( ! AspectFile.good() ) {
    cerr << "## GenerateLCP: Cannot open Aspect file " << AspectFileName
         << " ##" << endl;
    return false;
  }
  else
    SkipASCIIHeader( AspectFile, ElevHeader, AspectFileName, ElevFileName );

  fstream FuelFile( FuelFileName.c_str() );
  if( ! FuelFile.good() ) {
    cerr << "## GenerateLCP: Cannot open Fuel file " << FuelFileName
         << " ##" << endl;
    return false;
  }
  else
    SkipASCIIHeader( FuelFile, ElevHeader, FuelFileName, ElevFileName );

  fstream CoverFile( CoverFileName.c_str() );
  if( ! CoverFile.good() ) {
    cerr << "## GenerateLCP: Cannot open Cover file " << CoverFileName
         << " ##" << endl;
    return false;
  }
  else
    SkipASCIIHeader( CoverFile, ElevHeader, CoverFileName, ElevFileName );

  //Initialize new header.
  HeadData NewHeader;
  memset( &NewHeader, 0x0, sizeof(HeadData) );
  NewHeader.HUnits = NewHeader.BUnits = NewHeader.PUnits = NewHeader.DUnits =
  NewHeader.WOptions = 0;
  NewHeader.CrownFuels = 0;       // 1 have crown fuels, 0 if not
  NewHeader.GroundFuels = 0;      // 1 have ground fuels, 0 if not
  //-Set latitude degrees (negative for southern hemisphere).
  NewHeader.Latitude = Latitude;
  NewHeader.GridUnits = 0;        //DISTANCE UNITS: 0 meters, 1 feet
  //-FUEL OPTIONS
  // 0 no custom models AND no conversion file
  // 1 custom models BUT no conversion file
  // 2 no cust models BUT conversion file
  // 3 cust models AND conversion file needed
  NewHeader.FOptions = 0;
  NewHeader.EUnits = 0;       //ELEVATION UNITS: 0 metric, 1 english elevation
  NewHeader.SUnits = 0;       //SLOPE UNITS: 0 degrees, 1 percent
  //-ASPECT UNITS: 0 grass cats, 1 grass deg, 2 azimuth deg
  NewHeader.AUnits = 2;
  NewHeader.CUnits = 1;       //COVER UNITS: 0 cats(0-4), 1 percentages

  double HeightMult, BaseMult, DensityMult, DuffMult;
  fstream TreeHeightFile, CrownBaseHeightFile, CrownBulkDensityFile,
          DuffLoadingFile, WoodyFile;

  if( TreeHeightFileName.length() ) {
    //1=metric (m), 2=english (ft), 3=metric*10, 4= english*10
    NewHeader.HUnits = 1;
    HeightMult = 10.0;      // 10.0 if native units
    TreeHeightFile.open( TreeHeightFileName.c_str() );
    if( ! TreeHeightFile.good() )
      cerr << "## GenerateLCP: Cannot open Tree Height file "
           << TreeHeightFileName << " ##" << endl;
    else
      SkipASCIIHeader( TreeHeightFile, ElevHeader, TreeHeightFileName,
                       ElevFileName );

    NewHeader.CrownFuels = 1;
  }
  // 1 if metric (m), 2 if english (ft), 3 if metric*10, 4 if english*10
  else NewHeader.HUnits = 0;

  if( CrownBaseHeightFileName.length() ) {
    NewHeader.BUnits = 1;   // 1 if metric (m), 2 if english (ft)
    BaseMult = 10.0;        // 10.0 if native units
    
    CrownBaseHeightFile.open( CrownBaseHeightFileName.c_str() );
    if( ! CrownBaseHeightFile.good() )
      cerr << "## GenerateLCP: Cannot open Crown Base Height file "
           << CrownBaseHeightFileName << " ##" << endl;
    else
      SkipASCIIHeader( CrownBaseHeightFile, ElevHeader,
                       CrownBaseHeightFileName, ElevFileName );

    NewHeader.CrownFuels = 1;
  }
  else NewHeader.BUnits = 0;

  if( CrownBulkDensityFileName.length() ) {
    // 1 if metric (m), 2 if english (ft), 3 if m*100, 4 if eng*1000
    NewHeader.PUnits = 1;
    // 100.0 if in native units metric, 1000.0 if native units english
    DensityMult = 100.0;

    CrownBulkDensityFile.open( CrownBulkDensityFileName.c_str() );
    if( ! CrownBulkDensityFile.good() )
      cerr << "## GenerateLCP: Cannot open Crown Bulk Density file "
           << CrownBulkDensityFileName << " ##" << endl;
    else
      SkipASCIIHeader( CrownBulkDensityFile, ElevHeader,
                       CrownBulkDensityFileName, ElevFileName );

    NewHeader.CrownFuels = 1;
  }
  else NewHeader.PUnits = 0;  // 1 if metric (m), 2 if english (ft)

  if( DuffLoadingFileName.length() ) {
    NewHeader.DUnits = 2;
    DuffMult = 10.0;          // 10xMg/ha

    DuffLoadingFile.open( DuffLoadingFileName.c_str() );
    if( ! DuffLoadingFile.good() )
      cerr << "## GenerateLCP: Cannot open Duff Loading file "
           << DuffLoadingFileName << " ##" << endl;
    else
      SkipASCIIHeader( DuffLoadingFile, ElevHeader,
                       DuffLoadingFileName, ElevFileName );

    NewHeader.GroundFuels = 1;
  }
  else NewHeader.DUnits = 0;

  if( WoodyFileName.length() ) {
    NewHeader.WOptions = 1;

    WoodyFile.open( WoodyFileName.c_str() );
    if( ! WoodyFile.good() )
      cerr << "## GenerateLCP: Cannot open Woody file "
           << WoodyFileName << " ##" << endl;
    else
      SkipASCIIHeader( WoodyFile, ElevHeader,
                       WoodyFileName, ElevFileName );

    NewHeader.GroundFuels = 1;
  }
  else NewHeader.WOptions = 0;

  //Open output file.
  fstream OutFile( LandFileName.c_str(), ios::out|ios::binary );
  if( ! OutFile.good() ) {
    cerr << "## GenerateLCP: Cannot open Landscape file " << LandFileName
         << " for writing ##" << endl;
    return false;
  }


  memset( NewHeader.ElevFile, 0x0, sizeof(NewHeader.ElevFile) );
  memset( NewHeader.SlopeFile, 0x0, sizeof(NewHeader.SlopeFile) );
  memset( NewHeader.AspectFile, 0x0, sizeof(NewHeader.AspectFile) );
  memset( NewHeader.FuelFile, 0x0, sizeof(NewHeader.FuelFile) );
  memset( NewHeader.CoverFile, 0x0, sizeof(NewHeader.CoverFile) );
  memset( NewHeader.HeightFile, 0x0, sizeof(NewHeader.HeightFile) );
  memset( NewHeader.BaseFile, 0x0, sizeof(NewHeader.BaseFile) );
  memset( NewHeader.DensityFile, 0x0, sizeof(NewHeader.DensityFile) );
  memset( NewHeader.DuffFile, 0x0, sizeof(NewHeader.DuffFile) );
  memset( NewHeader.WoodyFile, 0x0, sizeof(NewHeader.WoodyFile) );

  NewHeader.XResol = ElevHeader.CellSize;
  NewHeader.YResol = ElevHeader.CellSize;
  NewHeader.numnorth = ElevHeader.NRows;
  NewHeader.numeast = ElevHeader.NCols;
  NewHeader.WestUtm = ElevHeader.XLLCorner;
  NewHeader.SouthUtm = ElevHeader.YLLCorner;
  NewHeader.EastUtm = ElevHeader.XLLCorner +
                      ( (double) ElevHeader.NCols*ElevHeader.CellSize );
  NewHeader.NorthUtm = ElevHeader.YLLCorner +
                       ( (double) ElevHeader.NRows*ElevHeader.CellSize );
  double ipart;
  modf( NewHeader.SouthUtm/1000.0, &ipart );
  NewHeader.lonorth = ipart * 1000.0;
  NewHeader.lonorth = NewHeader.SouthUtm - NewHeader.lonorth;
  modf( NewHeader.WestUtm/1000.0, &ipart );
  NewHeader.loeast = ipart * 1000.0;
  NewHeader.loeast = NewHeader.WestUtm - NewHeader.loeast;
  NewHeader.hieast = NewHeader.numeast * NewHeader.XResol + NewHeader.loeast;
  NewHeader.hinorth = NewHeader.numnorth * NewHeader.YResol
                      + NewHeader.lonorth;

  WriteHeader( OutFile, NewHeader );

  long   i, j;
  double maxval[10], minval[10];
  long   NumAllCats[10];
  long   AllCats[10][100];

  for( i = 0; i < 10; i++ ) {
    maxval[i] = -1e100;
    minval[i] = 1e100;
    NumAllCats[i] = 0;
    memset( &AllCats[i], 0x0, 100*sizeof(long) );
  }


  double elev, slope, aspect, fuel, cover, height, base, dense, duff, woody;
  CellData NewCell;
  CrownData NewCrown;
  GroundData NewGround;

  for( j = 1; j < NewHeader.numnorth + 1; j++ ) {
    for( i = 0; i < NewHeader.numeast; i++ ) {
      ElevFile >> elev;
      if( elev < 0.0 ) elev = -9999;
      NewCell.e = (short) elev;

      SlopeFile >> slope;
      if( slope < 0.0 ) slope = -9999;
      NewCell.s = (short) slope;

      AspectFile >> aspect;
      if( aspect < 0.0 ) aspect = -9999;
      NewCell.a = (short) aspect;

      FuelFile >> fuel;
      if( fuel < 0 ) fuel = 99;
      NewCell.f = (short) fuel;

      CoverFile >> cover;
      if( cover < 0.0 ) cover = 0.0;
      NewCell.c = (short) cover;

      WriteCell( OutFile, NewCell );

      if( NewHeader.CrownFuels ) {
        if( NewHeader.HUnits > 0 ) {
          TreeHeightFile >> height;
          if( height < 0 ) height = 0;
          height *= HeightMult;
        }
        else height = -1;
        NewCrown.h = (short) height;
        if( NewHeader.BUnits > 0 ) {
          CrownBaseHeightFile >> base;
          if( base < 0 ) base = 0;
          base *= BaseMult;
        }
        else base = -1;
        NewCrown.b = (short) base;
        if( NewHeader.PUnits > 0 ) {
          CrownBulkDensityFile >> dense;
          if( dense < 0 ) dense = 0;
          dense *= DensityMult;
          if( dense > 32767.0 ) dense = 32767.0;
        }
        else dense = -1;
        NewCrown.p = (short) dense;
        WriteCrown( OutFile, NewCrown );
      }
      else {
        NewCrown.h = 0;
        NewCrown.b = 0;
        NewCrown.p = 0;
      }

      if( NewHeader.GroundFuels ) {
        if( NewHeader.DUnits > 0 ) {
          DuffLoadingFile >> duff;
          if( duff < 0.0 ) duff = 0.0;
          duff *= DuffMult;
          if( duff > 32767.0 ) duff = 32767.0;
          NewGround.d = (short) duff;
        }
        else NewGround.d = -1;

        if( NewHeader.WOptions > 0 ) {
          WoodyFile >> woody;
          if( woody < 0.0 ) woody = 0.0;
          NewGround.w = (short) woody;
        }
        else NewGround.w = -1;

        WriteGround( OutFile, NewGround );
      }
      else {
        NewGround.d = 0;
        NewGround.w = 0;
      }

      FillCats( &NewCell, &NewCrown, &NewGround, AllCats, NumAllCats,
                maxval, minval );
    }
  }

  SortCats( AllCats, NumAllCats, minval );
  FinalizeHeader( &NewHeader, AllCats, NumAllCats, maxval, minval,
                  ElevFileName, SlopeFileName, AspectFileName, FuelFileName,
                  CoverFileName, TreeHeightFileName,
                  CrownBaseHeightFileName, CrownBulkDensityFileName,
                  DuffLoadingFileName, WoodyFileName );

  if( NewHeader.CrownFuels ) {
    if( NewHeader.HUnits == 5 ) NewHeader.HUnits = 1;
    if( NewHeader.BUnits == 5 ) NewHeader.BUnits = 1;
    if( NewHeader.PUnits == 5 ) NewHeader.PUnits = 1;
  }

  OutFile.seekp( 0, ios_base::beg );

  WriteHeader( OutFile, NewHeader );

  OutFile.close();
  ElevFile.close();
  SlopeFile.close();
  AspectFile.close();
  FuelFile.close();
  CoverFile.close();

  if( NewHeader.CrownFuels > 20 ) {
    if( NewHeader.HUnits ) TreeHeightFile.close();
    if( NewHeader.BUnits ) CrownBaseHeightFile.close();
    if( NewHeader.PUnits ) CrownBulkDensityFile.close();
  }
  if( NewHeader.GroundFuels > 20 ) {
    if( NewHeader.DUnits ) DuffLoadingFile.close();
    if( NewHeader.WOptions ) WoodyFile.close();
  }

  return true;
} //GenerateLCP

/*=============================================================================
  SkipASCIIHeader
  Skips and compares the ASCII header in a Arc ASCII grid file.
*/
void   SkipASCIIHeader( fstream& ASCIIFile, ASCIIGridHeader& Header,
                        const string& FileName1, const string& FileName2 )
{ //SkipASCIIHeaders
  string keyword, value;

  long CellSizeDiv2 = Header.CellSize / 2;

  for( int i = 0; i < 6; ++i ) {
    ASCIIFile >> keyword >> value;
    if( keyword[0] == 'n' || keyword[0] == 'N' ) {
      long val = atoi( value.c_str() );
      if( keyword[1] == 'c' || keyword[1] == 'C' ) {
        if( val != Header.NCols )
          cerr << "## SkipASCIIHeaders: NCOLS differs in " << FileName1
               << " from " << FileName2 << " (" << val << " versus "
               << Header.NCols << ") ##" << endl;
      }
      else if( keyword[1] == 'r' || keyword[1] == 'R' ) {
        if( val != Header.NRows )
          cerr << "## SkipASCIIHeaders: NROWS differs in " << FileName1
               << " from " << FileName2 << " (" << val << " versus "
               << Header.NRows << ") ##" << endl;
      }
      else if( keyword[1] == 'o' || keyword[1] == 'O' ) {
        //Assume this is the 'NODATA_VALUE' line.
      }
      else
        cerr << "## SkipASCIIHeaders: Odd header line in " << FileName1
             << ": " << val << " ##" << endl;
    }
    else if( keyword[0] == 'x' || keyword[0] == 'X' ) {
      double val = atof( value.c_str() );
      if( val < Header.XLLCorner - CellSizeDiv2 ||
          val > Header.XLLCorner + CellSizeDiv2 )
        cerr << "## SkipASCIIHeaders: XLLCorner differs in " << FileName1
               << " from " << FileName2 << " (" << val << " versus "
               << Header.XLLCorner << ") ##" << endl;
    }
    else if( keyword[0] == 'y' || keyword[0] == 'Y' ) {
      double val = atof( value.c_str() );
      if( val < Header.YLLCorner - CellSizeDiv2 ||
          val > Header.YLLCorner + CellSizeDiv2 )
        cerr << "## SkipASCIIHeaders: YLLCorner differs in " << FileName1
               << " from " << FileName2 << " (" << val << " versus "
               << Header.YLLCorner << ") ##" << endl;
    }
    else if( keyword[0] == 'c' || keyword[0] == 'C' ) {
      double val = atof( value.c_str() );
      if( val < Header.CellSize - 1 || val > Header.CellSize + 1 )
        cerr << "## SkipASCIIHeaders: CELLSIZE differs in " << FileName1
               << " from " << FileName2 << " (" << val << " versus "
               << Header.CellSize << ") ##" << endl;
    }
  }
} //SkipASCIIHeaders

/*=============================================================================
  FillCats
*/
void   FillCats( CellData *cell, CrownData *cfuel, GroundData *gfuel,
                 long AllCats[][100], long *NumAllCats,
                 double *maxval, double *minval )
{ //FillCats
  long k, m;

  AllCats[0][NumAllCats[0]] = cell->e;
  AllCats[1][NumAllCats[1]] = cell->s;
  AllCats[2][NumAllCats[2]] = cell->a;
  AllCats[3][NumAllCats[3]] = cell->f;
  AllCats[4][NumAllCats[4]] = cell->c;

  if( cfuel->h >= 0 ) AllCats[5][NumAllCats[5]] = cfuel->h;
  if( cfuel->b >= 0 ) AllCats[6][NumAllCats[6]] = cfuel->b;
  if( cfuel->p >= 0 ) AllCats[7][NumAllCats[7]] = cfuel->p;

  if( gfuel->d >= 0 ) AllCats[8][NumAllCats[8]] = gfuel->d;
  if( gfuel->w >= 0 ) AllCats[9][NumAllCats[9]] = gfuel->w;

  for( m = 0; m < 10; m++ ) {
    if( maxval[m] < AllCats[m][NumAllCats[m]] )
      maxval[m] = AllCats[m][NumAllCats[m]];
    if( AllCats[m][NumAllCats[m]] != NODATA &&
        AllCats[m][NumAllCats[m]] >= 0 )
      if( minval[m] > AllCats[m][NumAllCats[m]] )
        minval[m] = AllCats[m][NumAllCats[m]];
  }

  for( m = 0; m < 10; m++ ) {
    if( NumAllCats[m] > 98 ) continue;
    for( k = 0; k < NumAllCats[m]; k++ ) {
      if( AllCats[m][NumAllCats[m]] == AllCats[m][k] ) break;
    }
    if( k == NumAllCats[m] ) NumAllCats[m]++;
  }
} //FillCats

/*=============================================================================
  SortCats
*/
void   SortCats( long AllCats[][100], long *NumAllCats, double *minval )
{ //SortCats
  long i, j, m;
  long SwapCats[101];

  for( m = 0; m < 10; m++ )
    if( NumAllCats[m] > 98 ) NumAllCats[m]=-1;

  for( m = 0; m < 10; m++ ) {
    if( NumAllCats[m] < 0 ) continue;
    memcpy( SwapCats, AllCats[m], 100*sizeof(long) );
    for( i = 0; i < NumAllCats[m] - 1; i++ ) {
      for( j = i + 1; j < NumAllCats[m]; j++ ) {
        if( SwapCats[j] < SwapCats[i] ) {
          SwapCats[100] = SwapCats[i];
          SwapCats[i] = SwapCats[j];
          SwapCats[j] = SwapCats[100];
        }
      }
    }

    AllCats[m][0] = 0;
    for( i = 0; i < NumAllCats[m]; i++ )
      AllCats[m][i+1] = SwapCats[i];
    minval[m] = AllCats[m][1];  // zero slot is for nodata value
    if( minval[m] < 0 ) minval[m] = 0;
  }
} //SortCats

/*=============================================================================
  FinalizeHeader
*/
void   FinalizeHeader( HeadData *NewHeader, long AllCats[][100],
                       long *NumAllCats, double *maxval, double *minval,
                       string& ElevFileName, string& SlopeFileName,
                       string& AspectFileName, string& FuelFileName,
                       string& CoverFileName, string& TreeHeightFileName,
                       string& CrownBaseHeightFileName,
                       string& CrownBulkDensityFileName,
                       string& DuffLoadingFileName, string& WoodyFileName )
{ //FinalizeHeader
  char Description[512] = "";

  NewHeader->CrownFuels += 20;
  NewHeader->GroundFuels += 20;
  NewHeader->numelev = NumAllCats[0];
  NewHeader->numslope = NumAllCats[1];
  NewHeader->numaspect = NumAllCats[2];
  NewHeader->numfuel = NumAllCats[3];
  NewHeader->numcover = NumAllCats[4];
  NewHeader->numheight = NumAllCats[5];
  NewHeader->numbase = NumAllCats[6];
  NewHeader->numdensity = NumAllCats[7];
  NewHeader->numduff = NumAllCats[8];
  NewHeader->numwoody = NumAllCats[9];
  memcpy( NewHeader->elevs, AllCats[0], 100*sizeof(long) );
  memcpy( NewHeader->slopes, AllCats[1], 100*sizeof(long) );
  memcpy( NewHeader->aspects, AllCats[2], 100*sizeof(long) );
  memcpy( NewHeader->fuels, AllCats[3], 100*sizeof(long) );
  memcpy( NewHeader->covers, AllCats[4], 100*sizeof(long) );
  memcpy( NewHeader->heights, AllCats[5], 100*sizeof(long) );
  memcpy( NewHeader->bases, AllCats[6], 100*sizeof(long) );
  memcpy( NewHeader->densities, AllCats[7], 100*sizeof(long) );
  memcpy( NewHeader->duffs, AllCats[8], 100*sizeof(long) );
  memcpy( NewHeader->woodies, AllCats[9], 100*sizeof(long) );
  NewHeader->loelev = minval[0];
  NewHeader->hielev = maxval[0];
  NewHeader->loslope = minval[1];
  NewHeader->hislope = maxval[1];
  NewHeader->loaspect = minval[2];
  NewHeader->hiaspect = maxval[2];
  NewHeader->lofuel = minval[3];
  NewHeader->hifuel = maxval[3];
  NewHeader->locover = minval[4];
  NewHeader->hicover = maxval[4];
  NewHeader->loheight = minval[5];
  NewHeader->hiheight = maxval[5];
  NewHeader->lobase = minval[6];
  NewHeader->hibase = maxval[6];
  NewHeader->lodensity = minval[7];
  NewHeader->hidensity = maxval[7];
  NewHeader->loduff = minval[8];
  NewHeader->hiduff = maxval[8];
  NewHeader->lowoody = minval[9];
  NewHeader->hiwoody = maxval[9];

  //Copy all filenames into the header.
  sprintf( NewHeader->ElevFile, "%s", ElevFileName.c_str() );
  sprintf( NewHeader->SlopeFile, "%s", SlopeFileName.c_str() );
  sprintf( NewHeader->AspectFile, "%s", AspectFileName.c_str() );
  sprintf( NewHeader->FuelFile, "%s", FuelFileName.c_str() );
  sprintf( NewHeader->CoverFile, "%s", CoverFileName.c_str() );
  if( NewHeader->CrownFuels > 20 ) {
    if( NewHeader->HUnits > 0 )
      sprintf( NewHeader->HeightFile, "%s", TreeHeightFileName.c_str() );
    if( NewHeader->BUnits > 0 )
      sprintf( NewHeader->BaseFile, "%s", CrownBaseHeightFileName.c_str() );
    if( NewHeader->PUnits > 0 )
      sprintf( NewHeader->DensityFile, "%s",
               CrownBulkDensityFileName.c_str() );
  }
  if( NewHeader->GroundFuels > 20 ) {
    if( NewHeader->DUnits > 0 )
      sprintf( NewHeader->DuffFile, "%s", DuffLoadingFileName.c_str() );
    if( NewHeader->WOptions > 0 )
      sprintf( NewHeader->WoodyFile, "%s", WoodyFileName.c_str() );
  }
} //FinalizeHeader

/*=============================================================================
  WriteHeader
  Note that 'long' variables are written with a hard-coded size of 4. This is
  to enforce that the output file remains the same size, even if produced on a
  64-bit machine (instead of 32-bit).
*/
void   WriteHeader( fstream& OutFile, HeadData& Header )
{ //WriteHeader
  int    i;
  OutFile.write( (char *)&Header.CrownFuels, 4 );
  OutFile.write( (char *)&Header.GroundFuels, 4 );
  OutFile.write( (char *)&Header.Latitude, 4 );
  OutFile.write( (char *)&Header.loeast, sizeof(double) );
  OutFile.write( (char *)&Header.hieast, sizeof(double) );
  OutFile.write( (char *)&Header.lonorth, sizeof(double) );
  OutFile.write( (char *)&Header.hinorth, sizeof(double) );
  OutFile.write( (char *)&Header.loelev, 4 );
  OutFile.write( (char *)&Header.hielev, 4 );
  OutFile.write( (char *)&Header.numelev, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.elevs[i], 4 );
  OutFile.write( (char *)&Header.loslope, 4 );
  OutFile.write( (char *)&Header.hislope, 4 );
  OutFile.write( (char *)&Header.numslope, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.slopes[i], 4 );
  OutFile.write( (char *)&Header.loaspect, 4 );
  OutFile.write( (char *)&Header.hiaspect, 4 );
  OutFile.write( (char *)&Header.numaspect, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.aspects[i], 4 );
  OutFile.write( (char *)&Header.lofuel, 4 );
  OutFile.write( (char *)&Header.hifuel, 4 );
  OutFile.write( (char *)&Header.numfuel, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.fuels[i], 4 );
  OutFile.write( (char *)&Header.locover, 4 );
  OutFile.write( (char *)&Header.hicover, 4 );
  OutFile.write( (char *)&Header.numcover, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.covers[i], 4 );
  OutFile.write( (char *)&Header.loheight, 4 );
  OutFile.write( (char *)&Header.hiheight, 4 );
  OutFile.write( (char *)&Header.numheight, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.heights[i], 4 );
  OutFile.write( (char *)&Header.lobase, 4 );
  OutFile.write( (char *)&Header.hibase, 4 );
  OutFile.write( (char *)&Header.numbase, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.bases[i], 4 );
  OutFile.write( (char *)&Header.lodensity, 4 );
  OutFile.write( (char *)&Header.hidensity, 4 );
  OutFile.write( (char *)&Header.numdensity, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.densities[i], 4 );
  OutFile.write( (char *)&Header.loduff, 4 );
  OutFile.write( (char *)&Header.hiduff, 4 );
  OutFile.write( (char *)&Header.numduff, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.duffs[i], 4 );
  OutFile.write( (char *)&Header.lowoody, 4 );
  OutFile.write( (char *)&Header.hiwoody, 4 );
  OutFile.write( (char *)&Header.numwoody, 4 );
  for( i = 0; i < 100; ++i ) OutFile.write( (char *)&Header.woodies[i], 4 );
  OutFile.write( (char *)&Header.numeast, 4 );
  OutFile.write( (char *)&Header.numnorth, 4 );
  OutFile.write( (char *)&Header.EastUtm, sizeof(double) );
  OutFile.write( (char *)&Header.WestUtm, sizeof(double) );
  OutFile.write( (char *)&Header.NorthUtm, sizeof(double) );
  OutFile.write( (char *)&Header.SouthUtm, sizeof(double) );
  OutFile.write( (char *)&Header.GridUnits, 4 );
  OutFile.write( (char *)&Header.XResol, sizeof(double) );
  OutFile.write( (char *)&Header.YResol, sizeof(double) );
  OutFile.write( (char *)&Header.EUnits, sizeof(short) );
  OutFile.write( (char *)&Header.SUnits, sizeof(short) );
  OutFile.write( (char *)&Header.AUnits, sizeof(short) );
  OutFile.write( (char *)&Header.FOptions, sizeof(short) );
  OutFile.write( (char *)&Header.CUnits, sizeof(short) );
  OutFile.write( (char *)&Header.HUnits, sizeof(short) );
  OutFile.write( (char *)&Header.BUnits, sizeof(short) );
  OutFile.write( (char *)&Header.PUnits, sizeof(short) );
  OutFile.write( (char *)&Header.DUnits, sizeof(short) );
  OutFile.write( (char *)&Header.WOptions, sizeof(short) );
  OutFile.write( (char *)&Header.ElevFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.SlopeFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.AspectFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.FuelFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.CoverFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.HeightFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.BaseFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.DensityFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.DuffFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.WoodyFile, sizeof(char)*256 );
  OutFile.write( (char *)&Header.Description, sizeof(char)*512 );
} //WriteHeader

/*=============================================================================
  WriteCell
*/
void   WriteCell( fstream& OutFile, CellData& Cell )
{ //WriteCell
  OutFile.write( (char *)&Cell.e, sizeof(short) );
  OutFile.write( (char *)&Cell.s, sizeof(short) );
  OutFile.write( (char *)&Cell.a, sizeof(short) );
  OutFile.write( (char *)&Cell.f, sizeof(short) );
  OutFile.write( (char *)&Cell.c, sizeof(short) );
} //WriteCell

/*=============================================================================
  WriteCrown
*/
void   WriteCrown( fstream& OutFile, CrownData& Crown )
{ //WriteCrown
  OutFile.write( (char *)&Crown.h, sizeof(short) );
  OutFile.write( (char *)&Crown.b, sizeof(short) );
  OutFile.write( (char *)&Crown.p, sizeof(short) );
} //WriteCrown

/*=============================================================================
  WriteGround
*/
void   WriteGround( fstream& OutFile, GroundData& Ground )
{ //WriteGround
  OutFile.write( (char *)&Ground.d, sizeof(short) );
  OutFile.write( (char *)&Ground.w, sizeof(short) );
} //WriteGround
